Computing waveforms for spinning compact binaries in quasi-eccentric orbits 
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Several scenarios have been proposed in which the orbits of binary black holes enter the band of 
a gravitational wave detector with significant eccentricity. To avoid missing these signals or biasing 
parameter estimation it is important that we consider waveform models that account for eccentricity. 
The ingredients needed to compute post-Newtonian (PN) waveforms produced by spinning black 
holes inspiralling on quasi-eccentric orbits have been available for almost two decades at 2 PN 
order, and this work has recently been extended to 2.5 PN order. However, the computational 
cost of directly implementing these waveforms is high, requiring many steps per orbit to evolve the 
system of coupled differential equations. Here we employ the standard techniques of a separation of 
timescales and a generalized Keplerian parameterization of the orbits to produce efficient waveforms 
describing spinning black hole binaries with arbitrary masses and spins on quasi-eccentric orbits to 
1.5 PN order. We separate the fast orbital timescale from the slow spin-orbit precession timescale by 
solving for the orbital motion in a non-interial frame of reference that follows the orbital precession. 
We outline a scheme for extending our approach to higher post-Newtonian order. 



I. INTRODUCTION 

Standard scenarios for the formation of black hole bi- 
naries predict that the orbits will have circularized [T] by 
the time the system reaches the late inspiral phase probed 
by ground or space based gravitational wave detectors. 
There are however, alternative scenarios that produce 
systems with significant residual eccentricity for stellar 
mass systems 0, [1]; intermediate mass systems 0, 
and supermassive systems [6l-[ll|. Neglecting the effects 
of eccentricity on the waveforms will hurt detection \\3t - 
[T^ l and parameter estimation [TBI. 

In most scenarios the individual black holes that make 
up the binary will be spinning, so what is needed are 
waveform templates that include the effects of spin and 
eccentricity. The equations of motion describing such 
systems, along with expressions for the instantaneous 
waveforms and energy and angular momentum fluxes 
were computed to 2nd post-Newtonian (PN) order (order 
[v/cY in relative velocity of the system] by Kidder, Wise- 
man and Will in the early 1990's 0, Il7|: and more re- 
cently the calculation was carried to 2.5 PN order [T8l . [l9| . 
By evolving these expressions for the coupled ordinary 
differential equations describing the orbital motion, spin 
precession and energy and angular momentum decay, it is 
possible to generate waveforms that can be used for data 
analysis and parameter estimation. The drawback of this 
direct approach is the high computational cost associated 
with accurately capturing the rapid orbital motion along 
with the smaller and slower effects of periastron preces- 
sion, spin precession and orbital decay. A numerical im- 
plementation also loses control of the post-Newtonian ex- 
pansion, introducing higher order effects that can lead to 
(potentially spurious) chaotic behavior [20l - l22| . 

Our goal here is to develop an efficient approach for 
producing waveforms for spinning eccentric binaries with 
arbitrary spins and eccentricities for use in gravitational 
wave astronomy [2^ . Given the long history of the prob- 
lem, it is surprising that such waveforms do not already 
exist in the literature. Waveforms are available in the 



case of non-spinning eccentric binaries where analytic ex- 
pressions for the orbital motion and phase evolution are 
available at 3 PN and 3.5 PN order, respectively [23 - [30| . 
These can be combined with 2 PN accurate expression for 
the spinless waveform amplitudes [3l| to produce ready 
to use gravitational waveforms. Efficient waveforms are 
also known for circular binaries with spin out to 1.5 PN 
order for precessing systems, and out to 2 PN order for 
non-precessing systems [32|. For the more general case 
of eccentric binaries with spin, analytic expressions for 
the orbital motion and phase evolution have been com- 
puted to 1.5 PN order for equal mass systems and sys- 
tems with one spinning body [s^. These can be com- 
bined with the 1.5 PN order expression for the waveform 
amplitudes [s^] to produce ready to use waveforms for 
this sub-class of systems (however we caution that the 
expressions in Ref. (33j include some spurious spin-orbit 
effects that fail to vanish in the zero smn limit, and the 
waveform amplitudes quoted in Ref. [34| contain several 
errors. We discuss and correct both in what follows). 

Here we provide, for the first time, efficient waveforms 
at 1.5 PN order that describe eccentric binaries with 
arbitrary (but moderate) mass ratio and general spins. 
As is the case with the efficient waveform models de- 
scribed above, our derivation relies on the separation of 
timescales available in the problem that was first de- 
scribed by Damour [3^1 • Our work extends the 1 PN 
accurate Post-Keplerian parameterization of the orbits 
introduced by Damour & Deruelle [s^ to include the 
leading order, 1.5 PN spin-orbit effects, for both the ra- 
dial and angular motion (the radial component had pre- 
viously been computed at this order by Wex ^3)- The 
key new element in our derivation, which allows us to 
analytically compute both the radial and angular com- 
ponents of the post-Keplerian description of the orbital 
motion, is that we re-cast the equations of motion in 
a non-inertial frame that follows the spin-orbit induced 
precession of the orbit. This allows us to maintain the 
separation between the fast orbital motion (which can 
be solve analytically), and the much slower spin-orbit 
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precession (which must be solve numericaUy for general 
masses and spins). We go on to compute the 1.5 PN 
order spin effects on the secular decay of the eccentric- 
ity and semi-major axis, building on the 1 PN accurate 
treatment of Junker and Schafer [38| , and revisiting the 
earlier calculation of Gergely, Perjes and Vasuth [s^ , who 
used a different spin-supplementary condition. Our ex- 
pressions for the orbital decay match - after converting 
from the eccentricity and semi-major axis to energy and 
angular momentum - with the gauge invariant expres- 
sions for the orbital decay found earlier by Rieth and 
Schafer [4Cf|. We also revisit the calculation of the wave- 
form amplitudes 
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and correct several errors in the 
published expressions. Combining all these elements we 
arrive at ready to use, efficient, 1.5 PN order accurate 
waveforms for general, spinning, eccentric binaries. We 
conclude with some thoughts about continuing the cal- 
culation to higher post-Newtonian order. 



II. EQUATIONS OF MOTION 
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There are five constants of the motion: the magnitude 
of the angular momentum L, the individual spin magni- 
tudes 51, 5*2, the quantity L • Scfi, and the energy 
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The angular momentum has post-Newtonian and spin 
corrections: L = Ljv + Lpjv + Lso + • ■ • , which to 1.5 
PN order are given by 
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We begin by deriving a semi-analytic solution to the 
dissipationless equations of motion at 1.5 PN order. This 
we accomplish in two steps - first we find an analytic so- 
lution to the equations of motion in a non-inertial frame 
that precesses with the orbital plane. Next we calcu- 
late the time dependent rotation between the inertial and 
precessing frame using a fast-slow decomposition of the 
spin-orbit precession equations. 

The equations of motion at 1 PN order were solved 
by Damour & Deruelle [SG^ using a generalized Keple- 
rian parameterization of the orbits. In what follows we 
will focus on the 1.5 PN order corrections, which can 
be added to the earlier result to give the complete solu- 
tion at this order. The 1.5 PN equations of motion are 
most readily solved in the gauge defined by the Pryce- 
Newton-Wigner [41I, |42| (PNW) spin-supplementary con- 
dition. Suppressing 1 PN terms, the relative separation 
of the binary system r, the individual spins 81,82 and 
the orbital an gula r momentum L evolve according to the 
equations [33ll37| (in units where G ^ c~ 1) 
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Here M = mi+m2 is the total mass, fi = mim2/M is the 
reduced mass, L^v = /irxv is the Newtonian contribution 
to the orbital angular momentum and 
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with r/ = fi/AI. The momentum and velocity can be 
decomposed into radial and angular contributions: 
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where n = r/r, m — Ljv x n, u,. — r and p± — L/r. 
Combining the radial component of ([l]), p,- = fir, with 
© yields 



2E 



2M 

r 



2L- 8, 



(12) 



Squaring ([T]) and dropping higher order terms yields 
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Solving the above set of equations is complicated by 
the precession of the orbital plane, which causes v± to 
appear as a mix of azimuthal and longitudinal motion. 
The equations are more readily solved by transforming 
to a non-inertial frame that follows the precession of the 
orbital plane. The precessing frame is defined by the 
condition 
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The velocity in the precessing frame is given by 
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from which it follows that the orbital plane remains fixed 
in the precessing frame, L^v • fprc = 0, and the radial mo- 
tion is unchanged, Tpic = r- Introducing the azimuthal 
coordinate (j) in the orbital plane we have Vp^Q — 
and 
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The equations for f and can be solved by introduc- 
ing the generalized Keplerian parameterization of the or- 
bits [Hi: 
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where u is the eccentric anomaly, n = 27r/ is the mean 
motion with orbital frequency /, fc is the fractional peri- 
astron advance per orbit, a is the semi-major axis, and 
the regular Keplerian eccentricity has split into the triad 
of time et, radial e^, and angular e^, eccentricities . Here 
we have used a coordinate system defined by the Newto- 
nian angular momentum Lat and the line of sight vector 
to the source N: 
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Following Ref. [36] and restoring the 1 PN contributions, 
we find 
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For completeness we have included the 1.5 PN correction 
to the perihelion precession, fc, even though it is formally 
a 2.5 PN order term. 

The next step is to solve the spin-orbit precession equa- 
tions ([2]) to establish the time-dependent transformation 
between the inertial and precessing frames of reference. 
We begin by writing L^v = L -I- where L denotes the 
slowly changing, orbit averaged angular momentum, and 
(5L is a small periodic correction that varies on the orbital 
timescale. Adopting a similar decomposition for the two 
spins we find 
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where 



g{t) = 



e sm u 



(l-e2)3/2 (l-e2) 



1 



(32) 

In the above expression it is understood that we are using 
the Newtonian limit for u, cj), e and a. Note that the func- 
tion g{t) is periodic with period T — 1//. The (5L term 
causes a periodic variation in the observed waveforms 
of the same order as the 1.5 PN amplitude corrections 
discussed below. The slowly varying orbit averaged ex- 
pressions for the spins and angular momentum are found 
by numerically integrating the coupled set of differential 
equations 
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By solving the fast varying contribution to the preces- 
sion equations analytically we have reduced the compu- 
tational cost by a factor of (Af/a)'^/2 relative to solving 
the full equations. This completes our solution of the 
dissipationless motion. 

It is interesting to compare our solution to other ex- 
pressions in the literature. Our expressions for the 
quantities that enter the radial motion, a, 6^,64 and n, 
agree with those found by Wex [131, Konigsdorffer and 
Gopakumar [3l| and Keresztes, Mikoczi and Gergely ,43] . 
Konigsdorffer and Gopakumar [sl] provide expressions 
for higher PN order corrections that do not involve spin, 
while Keresztes, Mikoczi and Gergely [43] include the 
leading order spin-spin corrections. For the angular mo- 
tion Wex assumed that the spin precession could be ne- 
glected, leading to incorrect expressions for k and e^. 
Konigsdorffer and Gopakumar included the effects of spin 
precession, but their analysis was limited to the special 
case of simple precession, where either one spin vanishes 
or the bodies have equal mass. Our expressions for k and 
60 agree with theirs in the simple precession limit, save 
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for some spurious terms in their expressions that fail to 
vanish when the spins are set equal to zero. Keresztes, 
Mikoczi and Gergely did not provide a solution for the an- 
gular motion. More recently, Tessmer [3] has provided a 
solution for general spin orientations in the circular limit, 
but in a form that makes it difficult to compare to our 
solutions. 



III. DISSIPATION 

Dissipational effects first enter the equations of motion 
at 2.5 PN order. Rather than directly integrating these 
equations, we adiabatically evolve the system by incre- 
menting the energy and angular momentum according 
to the fiux equations. Because the dissipation occurs on 
a much longer timescale than the orbital or precession 
motion, we begin by orbit averagin g th e instantaneous 
flux equations. Junker and Schafer [38( 1 carried out the 
calculation to 1 PN order, and we now extend their cal- 
culation to include the spin-orbit effects at 1.5 PN order. 
The instantaneous expressions for the 1.5 PN fluxes were 
computed by Kidder [17[ using the covariant spin supple- 
mentary condition, and by Zeng and Will [45j using the 
PNW spin supplementary condition. We use the latter 
expressions to be consistent with the choice we made for 
the dissipationless equations of motion. 

The orbit- averaged flux at 1.5 PN order has two con- 
tributions: one from averaging the PN flux over a 1.5 
PN order orbit, (-^0)1.5; the another from averaging the 
1.5 PN flux over a PN order orbit, (^1.5)0. For the 
energy these are: 
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Here we have introduced the spin combinations 
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which are related to Scs by 
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The contributions to the decay of the angular momentum 
are: 
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(332 + 15726^ + 4356-*) 



(42) 



where S|| = L(L • S), = S — Sy, and similarly for 
Z|| and Z±. The terms parallel to the orbital angular 
momentum contribute to the decay of the orbit (since 
i = L • L). The orthogonal terms introduce additional, 
higher order precessional effects. Rieth and Schafer [i^l 
found expressions for the spin-orbit corrections to the 
energy and angular momentum loss using a clever trick 
that allowed them to solve for the orbit averaged fluxes 
without having to solve the equations of motion. They 
quote their expressions in terms of E and L, but we can 
invert (j25p to flnd expressions for E and L in terms of a 
and Cr'. 
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and after substituting these expressions into their equa- 
tion (61) and converting the result into an energy flux by 
dividing through by — /z, we get perfect agreement with 
the sum of the terms in our equations and ([57)1 . Sim- 
ilarly, after correcting their equation (70) by restoring a 
missing factor of 1/L^ in the second term, we get perfect 
agreement with our expression for the rate of change of 
the angular momentum from equations (HIT) and (|42p . 
It is worth noting that the expressions for the orbit- 
averaged energy and angular momentum loss quoted by 
Rieth and Schafer are gauge invariant (i.e. they do 
not depend on the choice of SSC), while our expressions 
are not gauge invariant since they involve non-invariant 
quantities such as a and 6^. 

The adiabatic decay of the orbits is found by applying 
the chain rule to our expressions (P5|) for the semi-major 
axis and radial eccentricity: 
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and similarly for a. The eccentricity evolution provides 
a good cross check for our calculations as the collection 
of terms in square brackets have to cancel to order 62 to 
avoid unphysical behavior in the circular limit. 
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Putting everything together and restoring the PN 
and 1 PN terms we have 
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IV. WAVEFORMS 
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Gergely, Perjes and Vasiith [39'| have performed a simi- 
lar calculation using the covariant SSC, and while it is 
generally impossible to compare results that use differ- 
ent choices for the SSC, we note that the terms involving 
L • S are in agreement, while those involving L • Z are not. 
The explantion is that the L • S terms in the equations 
of motion are the same in the two SSC frames, while the 
L • Z terms differ. 



Our solution to the equations of motion is completed 
by using (|43p to recast the other Keplerian parameters 



With the orbital motion in hand, the final step is to 
derive expressions for the polarization states of the grav- 
itational waves. These can be written as 



h+ = ^{PiPj - qiqj)h'4^ 



(52) 



where hl^r^ are the transverse-traceless components of the 
metric perturbation and p and q are unit vectors orthog- 
onal to the line of sight vector N defined by 



p = Ljv X N/ILat X N| , 
q = N X p . 



(53) 



Expressions for hlj^^ out to the requisite order have 
been computed by Kidder [l3] using the covariant spin- 
supplementary condition. To our knowledge, a similar 
calculation has not been done using the PNW spin sup- 
plementary condition, but transforming Kidder's expres- 
sions to the PNW frame reveals no change to 1.5 PN 
order, which is all we need for our current purposes (the 
first changes appear at 2.5 PN order). 

Majar and Vasiith [s^] have provided formal expres- 
sions for the instantaneous polarization states using a 
corotating coordinate system. They first define an "in- 
variant" source frame using the line of sight vector N 
and the total angular momentum vector J = L + S, 
which they then relate to a corotating system attached 
to the orbital angular momentum Ljv and orbital sepa- 
ration vector r. When dissipation is included, J slowly 
evolves and it is necessary to introduce a new invari- 
ant coordinate system. A convenient choice might be a 
frame related to the detectors, such as the Barycenter or 
Geocenter frames, which are approximately invariant on 
the observational timescale of a space based or ground 
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based detector. Thus our description of the waveforms 
involves three reference frames: the detection frame; the 
source frame; and the orbital frame. The rotation be- 
tween the source frame and orbital frame occurs on the 
orbital timescale, while the rotation between the source 
frame and the detection frame occurs on the more sedate 
dissipation timescale. 

The source frame is described by the triad {i,j,k}, 
where k = J, i=jxJ, j=Jx N/sin7, and cos 7 = 
J • N. The comoving orbital frame is described by the 
triad {x, y,z}, where x = f , z = Ljy and y = Ljy x f. 
The corotating system is related to the source frame by 
three time dependent Euler angles {$, is, ^'}, where is 
the orbital phase in the source frame, $ is the precession 
angle of the orbital plane, and is is the precession cone 
angle. Applying these rotations to N, we find it has the 
following components in the corotating frame: 

Nx = cos 4* cos $ sin 7 — sin (sin $ sin 7 cos is 
— cos 7 sin 15) 



and 



sin cos $ sin 7 — cos 5'(sin<I' sin 7003*5 



— cos 7 sin 15) 
Nz — sin 15 sin $ sin 7 + cos 7 cos is . 
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The components of p and q in the corotating frame are 
then 



Pk ^ 0, 
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qk = + ■ 
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The mapping between the detection frame and the source 
frame can be accomplished by evolving the components of 
the spins and angular momentum in the detection frame, 
then solving for the angles that appear in the corotating 
frame using the relations 



cos is 
cos 7 

cos $ 



J Ln 
J-N 
N X J 



(^J X Ljv^ 



IN X JIU X L 



(57) 
(58) 

(59) 



and 132 



cos is $ . 
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For completeness, and also to correct several errors in 
the original paper, we re-compute and quote the explicit 
expression for the polarization states out to 1.5 PN order. 
The 1.5 PN terms are separated into spin-orbit (SO) and 
other contributions. 
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(7 + 3r])M _^ (1 - 5r?)f ^ 



(3-8?7)M (l-5?7)?; 
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4r 



(65) 



hi 



M 
r 



PxQx 



6m 
M 





'4M „.2' 








[(^ 


2f2 




2f^ 


Ny^ 




r 




r 



3M 



2r' 



-/Vx - 2v±rNy^ {p^Qy + q^py) 



-2v']_ {rN^ + v±Ny) pyqy 



hi = 



3((l-3r?) 



+ 



+ 6P + 



7M2 



+4v_Lr 



6M 



r 



n: 



+ 



(19 - Qrjy'^M 



+ (3 - 9r;)t>^r 



2^2 (10 + 3r])v'^M 29M=^ 



(2 + 477) M 



+ 6r 

+ |^((l-3,?)(2 





5M „.2' 








+ 2r^ 


Nl + d,v^ 


-A— + 3f=^ 




r 




r 



Pxqx 
N^Ny + 12vlfN\ 



+ (1 - 37?)z;2 
7M 



(Pa;gj/ + qxPy) 



+ 3f" 



(4 - 6r;)M 



(3 - 9r/)t;^ 



lSO 



-^[(A.q)g^-(A.p)p^] , 



5m 



(1 - 277) r 



SPM 2t;2M 41M2 



2r r 



g^2 



2f^ 



4r 



+ 6rv^ 



'2M .2' 




M , 




r 




2r^ 

2r 





+ v± 



(10 + 77?)f2M _^ (2 + 77)t;^M _ (59 - 307?)M^ _ _ ^^^^2^2' 



r r 6r^ 

(25 + 26r])r^M (7 - 277)t;2M (26 - 3r])M'^ 



+v± 



dm 
M 



4r 



(1 - 27?) 



4r 



+ 37;1 



5M 
'2r 



-2r' 



Ar 4r 

r2 o„,3 J. Ar3 



3r2 



2f^ 



(1 - bri)v^f 



\ Pxqx 







Nl + vi_f 


r 



NlNy 



(32 - 97?) M2 



3r2 



+v 



6m 



(1 - 27?) 



(1 - 5r])v'^f^ 
5M 



(49 + 147?)r-2^f (11 - 6r])v'^M 

Ar Ar 
■(2 + 67?)M 



2r 



(14 + 67?)M 



+ (1 - 5r/)'6'2 



Ax + v± 



r 

29M 



+ (1 - 57?)^^ 



Ny \ {jp^qy + qxPy) 



2r 

(3 - 8?7)M 



NlNy - QvlfNxNl - 2vlNl 



2r 



(1 - 57?)v^ 



Ny \ Pyqy , 
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where 5m = m2 — mi, and A = M(S2/to2 — Si/mi). 

We have checked that the waveforms we get by sub- 
stituting in our solution for the orbital motion and spin- 
orbit precession are in perfect agreement with the Lang 
and Hughes waveforms [46j] in the circular limit. 



V. FUTURE WORK 

It would be desirable to extend our treatment to higher 
post-Newtonian order. The necessary building blocks are 
known to 2.5 PN order, and while there is no fundamen- 
tal barrier to going to higher order, there are some new 
effects to contend with. We will briefly describe some of 
the issues that crop up at 2 PN order. At this order the 
precession equations read 



3 
3 



5cff X Lat - — ((n • S2)(n x Si) 

-f(n.S2)(nx Si)), (71) 



from which it follows that L • Scg and the are no 
longer constant. The condition (iLjv/c?i|prc = demands 
that the precessing frame rotates with angular velocity 

3n 



dSi 




IT 








dt 




dliN 


1 


dt 





Scff 



((n-S2)(I 



'N 



Si) 



n.Si)(LAr.S2) 

(72) 

Note that the new terms in w have no affect on the map- 
ping between the velocity in the inertial and precessing 
frames. As before, the equations for f and (j) can be writ- 
ten as polynomials in 1/r, with new spin^ terms from the 
2PN order Hamiltonian. In contrast to what we found 
at 1.5 PN order, the spin-dependent coefficients in the 
polynomial expansion are no longer constant. The time 
dependence of the coefficients prevents us from finding 
an exact generalized Keplerian solution to the equations 
of motion in the precessing frame beyond 1.5 PN order. 
On the other hand, we know that the spin dependent co- 
efficients are approximately constant on the orbital time 



scale since the time dependence enters at 2 PN order. 
Thus we can find a solution for the Keplerian parameters 
by treating the spin dependent terms as constants, which 
are then updated adiabatically via the precession equa- 
tions, just as the total energy and angular momentum 
are updated adiabatically via the dissipation equations. 

As first noted by Damour and SchaferfJ?!, and later 
solved by Schafer and Wex [48], the generalized Keple- 
rian parameterization has to be further generalized to 
handle 2 PN terms in the equations of motion for non- 
spinning bodies. In particular, they found it necessary 
to introduce terms involving the quantity 



1 



v^2tan-' (±^) tan^ , (73) 
\1 - e^y 2 

which is closely related to the orbital angle (p. Some of 
the spin-spin terms that appear at 2 PN order will require 
addition terms that depend on v. In particular, there are 
terms in the equation of motion of the form (Si-n)(S2-n), 
which produce terms like (Si • p)(S2 • p) cos^ v in the gen- 
eralized Keplerian parameterization. In the adiabatic ap- 
proach, spin dependent coefficients such as Si • q, S2 • p 
and Si • S2 are treated as constants when solving the 
equations of motion. Expressions for the 2PN spin-spin 
corrections to the radial motion have already been de- 
rived using the PNW spin supplementary condition [i^ . 
However, given that the terms involving L • Scff are no 
longer constant when 2 PN effects are taken into account, 
there is little advantage to using the PNW spin supple- 
mentary condition beyond 1.5 PN order. Indeed, it may 
be wise to adopt the covariant spin supplementary con- 
dition since this is the gauge in which the higher order 
corrections to the instantaneous waveforms have already 
been derived. 
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